#!/usr/bin/perl
# 
#  Exercise 8.9 
# 
# Write a program that randomly mutates the amino acids in a protein but restricts the 
# possibilities to those that can occur due to a single mutation in the original codons, 
# as in Exercises 8.7 and 8.8.
#
# Answer to Exercise 8.9
#
# We'll use the answer to Exercise 8.8 for this.  Because that answer did not always
# change the amino acid (which often happens when a codon is mutated to another codon
# that codes for the same amino acid), the peptide may not change at any given step
#

use strict;
use warnings;
use BeginPerlBioinfo;

# Call srand (may not be necessary if you're using a recent version of Perl
srand(time|$$);

print "Give a peptide (string of amino acids): ";
my $peptide = <STDIN>;
chomp $peptide;

#
# Use the subroutine from Exercise to ensure this is a protein.
#

die "$peptide is not a protein sequence\n" unless isProtein($peptide);

my $number_of_mutations = 10;

for(my $i=0 ; $i < $number_of_mutations ; ++$i) {

	my $position = rand length($peptide);

	my $new_aa = random_mutate_aa(substr($peptide,$position,1));

	substr($peptide,$position,1,) = $new_aa;

	print "$peptide\n";
}

exit;


##################################################
# Subroutines
##################################################

sub random_mutate_aa {

	my($aa) = @_;

	#
	# First calculate the underlying codon(s)
	#

	my(%reverse_genetic_code) = (
	
	'A' => 'GCA GCC GCG GCT',		# Alanine
	'C' => 'TGC TGT',			# Cysteine
	'D' => 'GAC GAT',			# Aspartic Acid
	'E' => 'GAA GAG',			# Glutamic Acid
	'F' => 'TTC TTT',			# Phenylalanine
	'G' => 'GGA GGC GGG GGT',		# Glycine
	'H' => 'CAC CAT',			# Histidine
	'I' => 'ATA ATC ATT',			# Isoleucine
	'K' => 'AAA AAG',			# Lysine
	'L' => 'CTA CTC CTG CTT TTA TTG',	# Leucine
	'M' => 'ATG',				# Methionine
	'N' => 'AAC AAT',			# Asparagine
	'P' => 'CCA CCC CCG CCT',		# Proline
	'Q' => 'CAA CAG',			# Glutamine
	'R' => 'CGA CGC CGG CGT AGA AGG',	# Arginine
	'S' =>'TCA TCC TCG TCT AGC AGT',	# Serine
	'T' => 'ACA ACC ACG ACT',		# Threonine
	'V' => 'GTA GTC GTG GTT',		# Valine
	'W' => 'TGG',				# Tryptophan
	'Y' => 'TAC TAT',			# Tyrosine
	'_' => 'TAA TAG TGA',			# Stop
	);

	my @possible_codons = split(' ', $reverse_genetic_code{$aa});

	#
	# Next, get the list of possible amino acids
	#

	my @possible_amino_acids =  ();

	foreach my $codon (@possible_codons) {

		push(@possible_amino_acids, codon_mutated_and_translated($codon));
	}

	#
	# Finally, randomly select one and return it.
	# Note that the array @possible_amino_acids may include some amino acids
	# more than once.  How would you remove all duplicate amino acids from
	# this array?  (One suggestion: sort the array, then loop through it
	# looking for adjacent duplicate amino acids.  Another suggestion:
	# use the trick in the subroutine codon_mutated_and_translated that
	# avoids duplicate entries in a hash.)
	#

	return $possible_amino_acids[rand @possible_codons];
}

#
# Here are subroutines from Exercise 8.7 we need
#

sub codon_mutated_and_translated {
	
	my($codon) = @_;

	my %translated_codons = ();

	my @mutated_codons = mutate_codons($codon);

	# A debugging print statement
	# print "The mutated codons are @mutated_codons\n";

	foreach my $mutated_codon (@mutated_codons) {

		# A debugging print statement
		# print "XX $mutated_codon -> ",codon2aa($mutated_codon),"\n";

		# This hash statement has the effect of defining each amino acid only once, so
		# the return statement for the subroutine doesn't give an amino acid twice if
		# there are two mutations that result in that amino acid.

		$translated_codons{codon2aa($mutated_codon)}++;
	}
	
	return keys %translated_codons;
}

#
# A helper subroutine: given a codon, it returns all the codons that
# can result from a single base change -- 9 in all.
#

sub mutate_codons {

	my($codon) = @_;

	my @bases = ('A','C','G','T');

	# To collect the output mutations
	my @mutations = ();
	
	#
	# Note the loop within a loop.  The outer loop loops through the 3 positions in the
	# codon.  The inner loop loops over the possible bases.
	#
	for(my $i = 0; $i < 3 ; ++$i) {
		foreach my $base (@bases) {

			my $mutation = $codon;

			if ($base eq substr($codon,$i,1)) {
				next;
			}else{
				substr($mutation,$i,1) = $base;
			}
			push(@mutations, $mutation)
		}
	}
	return @mutations;
}

# From Exercise 8.10
# 
sub isProtein {

	my($sequence) = @_;

	# If there's a non-amino acid (case-insensitive) character, return "false"
	if($sequence =~ /[^ACDEFGHIKLMNPQRSTVWXYZ]/i) { return 0; }

	# Otherwise, return "true"
	return 1;
}
